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Thic paper concerns a numerical model of a regenerator running at very lew 
tempe rat ores* The model consists of the usual three equations for a 
corspreasible fluid with an additional equation for a "matrix temperature, " The 
main difficulty with the model is the very low Mach number (approximately 
l.E-3)* The divergence of the velocity is not small, the pressure divergence is 
small, and the pressure fluctuation in time is not small* An asymptotic 
expansion based on the "bounded derivative" method of Kreiss is used to give a 
"reduced" model which eliminates acoustic waves* The velocity is than 
determined by a two~poin£ boundary value problem which do£3 not contain e time 
derivative. The solution obtained from the reduced system ic compared v;ith the 
numerical solution of the original system* 
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1* The differential equations 

The oscillating flew tskeG place in a cylinder filled with email metal spheres [lj. The 
compressible flow equations are used to describe the flat; around the spheres* A "resistance" terra 
involving the resistance factor F is used to account for the porous acedia effects* The "matrix" of 
Spheres’ h’ac a heat capacity and neat can be transferred between the matrix and the gaa. This 
conductance is determined by the coefficient h. The temperature of the matrix is given by T m * The 
variables are listed in the appendix. The partial differential equations for these variables are 
the following. 
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The domain io 0 £ x jC L. We will discuss the boundary and initial conditions in the next section. 
Their determination is crucial to the success of the numerical approximation. 

To set up the numerical approximation we need to write these equations in dimensionless form. 
For this purpose we assume an ideal gas described by the following relatione. 


p« pc* « - K c - 2.03 E+3 J/Jtg*K y - 0.6 

R C T Y 


The scaling factors for the bade variables are the constants . (p, u, T). The time is scaled by the 
period (t) of the imposed oscillation in the mass flow. This is the some ae the period of the 

piston oscillation. The x coordinate is scaled by the length L. We use (p, u, T, T ra ) for the 
scaled variables (p/p, u/u, T/dT, T m / AT). 
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Typical values of the 

oe dimonslonleoQ parameters 

are: 



*1 * 

0.92 

T 2 “ 

14 

«1 - 

70 

a 2 - 709 

03 ■ 

64 
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5.82-4 

a 2 « 

8.7E-3 

c - 1.4E-5 
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on the 
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t « 

0.5 o 

U « 

0.12 e/s 

L *» 

0.057 ta 

R c » 2.03B+' 

®r - 

0.65 

Cv “ 

3. 1E+3 j/kg*K 

h « 

8.7E+3 W/o^K 

- 


15 K 


Dj* « 1.2E-4 


Clearly the term with the dominant coefficient is the one containing the factor e~*. The 
pressure cuet be very nearly constant in order to balance the terras in theee equations. In order to 
avoid high velocity sound waves wo suet carefully set the initial and boundary conditions. The 
"bounded derivative" method of Rreiso io used to set the initial conditions. We chose the initial 
conditions so that the first two derivatives of the basic variables with respect to time are 0 ( 1 ) 
rather than 0(e~*). The work of Kreisa indicates that we can expect the time derivatives to remain 
bounded during the time Integration of the partial differential equations. This implies that the 
fast sound waves do not appear, since their presence would require derivatives which are 0 (e~*). 

2. The bounded derivative expansion 


This method wes developed by Krei30 [4]» Gustafoeon has shewn that, in certain cases, the 
results of Kreins can be obtained by an asymptotic expansion [3). Guotafeson has given a good 
expository treatment of the method [2). If we assume an expansion of the basic variables iti the fora 


p (x, t ) « p 0 (x,t) + epj(x,t) + ••• 


then we obtain a new cystem of equations for the functions (p<, T^, u^, T^), i * 0 , 1 , ... by 
equating like powers of e in the usual way. Thus, we obtain from the velocity eq (7) 
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The first term in the pressure expansion is therefore independent of x. Using this fact, eq (1) can 
be differentiated with respect to x to obtain an equation for the velocity (hereafter we drop the 
subscript "o" and the **' k " modifier for the asymptotic expansion). 
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If the functions T and p are known and u is known at the boundary, then this equation can be solved 
for u(x) in the interior. Equation (9) implies that the right side of the pressure eq (1) is 
independent of x. Therefore eq (1) can be regarded as an ordinary differential equation for p(t). 


201 



Equation (2) can be regarded as a hyperbolic partial differential equation for T(t,x) and 
eq (A) is an ordinary differential equation for T n (t,x). Note that these equations do not require 
any boundary conditions for the pressure p or the temperature T m . If we regard the equation for the 
gas temperature as a hyperbolic equation, then we must specify T at an inflow boundary and leave it 
unspecified at an outflow boundary. He want to specify constant temperatures T(t,o) ** and 
T(t,L) “ at the boundary. However, when the flow reverses this would cause a discontinuity in 
the temperature at the boundary. Therefore we impose a fast exponential decay in time from the 
temperature at the time of reversal to the desired temperature Tl or 

The velocity boundary values are chosen so that the mass flow is sinusoidal, that is 


p(t,o)u(t,o) a C 0 ein(wt) 
p(t,L)u(t»L) - 0^ ein(wt+8) 


These conditions insure that mass is conserved over a cycle. 

3. The numerical scheme for the asymptotic expansion 

A finite difference meah (t n ,X|) is used where t n « nAt, and jtj « iL/N for 0<i<N# The opatlal 
derivatives are approximated by centered second order differences, for example 
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where la the approximation for T(t n ,X|). An implicit type of predictor corrector is used which 
is similar to a Crank-Nicolaon scheme. We will not write out the complete scheme for the full 
equations; instead we describe the scheme for the following simple equation 
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Given the values at the nth time level, U^, and an approximation for , then compute a 
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corrected approximat ion Uj from 
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This is a tridiagonal linear systeu frr the unknown vector Uj. A correction is made for the 
variables p,T, and T n in that order; then an updated value of the velocity is obtained by solving a 

n 

finite difference version of eq (9) for U^. If the heat transfer cofficient, h, is constant, then 
this is another tridiagonal linear system. If h - h(u) depends on the velocity, then the resulting 
nonlinear equation Is solved by a Newton Iteration. The boundary conditions are described in the 
previous section. Generally, from two to five iterations (or corrections) were used. On the first 

~o n 
iteration « u^. 
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A. A fully Implicit scheme for the original equations 

have a second computer code which is based on the original eqs (1) through (A) rather than 

2 

the asymptotic expansion. Since. M 0 ° e is small we expect this sytem to be very stiff. Therefore, 
we expect that an implicit schema is required. However, we use a CrJ (Crank-Nicoleon) scheme rather 
than a DDF (backward-dif ference-forraula) . We choose the initial conditions to avoid the fast moving 
acoustic waves; therefore we do not need the BDF scheme to damp out the fast waves. We can choose 
arbitrary initial values for T(x,o); however p(x,o) is required to be constant and u(x,o) caist 
satisfy (9). This is the bounded derivative principle of Kxeise. The initial conditions are chosen 
so that the first two time deriativea of the solution (p, T, u, T B ) at t ■ 0 Are bounded independent 

Of E* 

The finite difference scheme uses three-point centered approximations for the spatial 
derivatives except at outflow boundaries, where a one-sided first order approximation is used. A 
Crank-Nicoleon approximation is used in time. The difference scheme for each of the four variables 
is thus similar to (1) for the asymptotic approximation. However, the four equations are now 
coupled. Therefore, each time step requires the solution of a linear block-tridiagonal system with 
AxA blocks. We do not use a Newton iteration to deal with the nonlinearity* Instead, we use an 
iteration and evaluate the coefficients of the derivative terms at the previous iteration. 

The boundary approximation for the temperature T(x,t) is the same cs for the asymptotic 
equations. The boundary conditions for u(x,t) are specified at all times using the oaoe values as 
in eq (6). Thera is no boundary condition for the p equation; instead one-sided differences are . 
used to approximate derivatives at the boundary. No boundary condition is needed for the T ffi 
equation since this equation contains no spatial derivatives. 

5. Computational results 

In this section we give some results obtained from the two numerical methods described in the 
previous sections. Unless otherwise stated, the results were obtained using the bounded derivative 
(t.e. asymptotic expansion) model. Within SI, many units choices exist for most quantities! In all 
of these runs the heat transfer coefficient between the gas and matrix was given by the term 


h(m) * AO cxp(-i.6m) + 220( I~exp(~l . 6m) )s> 


where 


m ° p|u| . 


The exponential factor is Included to avoid a discontinuous derivative of h(|pu|) with respect to to 
u at u*=0. 

The formula used to bring the boundary temperature of the gas back to its constant inflow value 
when the velocity reverses is 


T(t) + y^rev ““ Tijexp^t rev -t^/t j 


Here T^ is the constant inflow temperature (15 K at the left boundary and 10 K at the right boundary 
for most of our rune) and T rev is the temperature of the gas at the tine when the velocity reverses 
direction. The time of reversal i3 t rev and x Is an input parameter, whose value was G.05P where P 
is the period of the macs flow oscillation. 

The mass flow at the boundary is given by 


pu 13 C sin(2Trft+0) 
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where f ie the frequency and 0 la the phage. At the left boundary 0*0, at the right values In the 
range -45°<0<O° vere used. 

All runs assumed an Ideal gas with equation of otate 
p « Rc pT. 

The specific heat of the gas at constant volume was 
C v - 3120. J/kg’K 


The thermal conductivity of the matrix was 0.3 W/m*K and of the gas was 0.13 W/m*K. The regenerator 
length was 0.0572 m. 


5.1 The temperature and pressure wave forma 


The temperature as a function of time at the left and right boundaries is given In figures l 
and 2. In this case, which we refer to as case I, the inflow temperatures at these boundaries were 
15 K and 10 K; the mass flow amplitude wa3 77.9 kg/rrs; the starting pressure was 3.3 Mpa; the heat 
capacity of the matrix was 2.64E+5; the hydraulic diameter was 1.21E-4; and the frequency was 5 IU. 
The time interval shewn Is 12Ct<13, that is, tne system has been run for GO cycles before these 
curves are drawn. The matrix temperature is shown as a dashed line; the gas temperature is a solid 
line. In this esse, the two temperatures were almost identical. The temperature at a p^lut 4/5 of 
the distance across the regenerator (x^O.G) is ohc*m In figure 3. 

The temperature as a function of position k across the regenerator at certain times in the 
cycle is shown in figure 4. The pressure oscillation is shown in figure 5. 



Figure l. Temperature (K) vs. time <n) at x -.0 for case I. 
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Figure 2, Temperature vs* time at x ° 1.0 for case I. 



Figure 3. Temperature vs. time at x =• 0.8 for case I. 
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Figure 7. Temperature vs. time at x = 1.0 with p m C m = 5.28E+6 J/kg'ra^. 


mow some separation of the matrix and gas temperature. A greater reparation is seen in figure 7. 
Here the heat capacity of the matrix has been increased by a factor of 20 over case I. 

Our definition of the ineffectiveness is given by the following integral which a taken over 
the outflow portion of the cycle, 


/ p |u|(Ta,tW U )dC / 7p|u|(T 0 -T L )dt. 


Here T Q and T L denote the fixed gas temperature and T(L,t) denotes the gas temperature o . outflow. 
The ineffectiveness as a function of hydraulic diameter, D^," is shown in figure 8 for the conditions 
of case I, In figure 9 the variation with matrix heat capacity, Pro C m » ** s s^ 0 *™ 4 

In table 1, va compare our computed ineffectiveness with that obtained by Daney using a model 
which assumes no pressure variation. Our results are higher, which is probably explained by our 
non-zero pressure swing during the cycle. The table shows that ineffectiveness increases with the 
pressure swing. 


5.3 The accuracy and consistency of the model 

In table 2 the value of the ineffectiveness is shown as a function of the numerical resolution. 

These results are for case I where the frequency is 5 Hz . Note that with At « 0.001 s and Ax - 0.1 
we have 10 mesh intervals across the regenerator and 200 time steps per cycle. The last run with At 
* 0.004 s did net yield valid results. The need for a small time step may be explained by the rapid 
change in the temperature at the boundary. When the flaw changes from outflow to inflow the- 
tenv-srature ie raised to the boundary value in about 0.01 seconds in this case. Therefore we might 
expect some trouble with the larger time step. 
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HYDRAULIC DWWETER, 

f fectiveoecs vs. hydraulic diameter, 




Tabic 1* CoEparison with Constant Pressure Model- 


Model 


Temperature Swing 
Ineffectiveness K 


Pressure Swing 
MPa 


constant pressure (Daney) 0.015 
phase ° 0.0 0.027 
phase • -14° 0.079 
phase - -29° 0.122 


10.0 - 20.50 2.9 - 3.3 
10.0 - 10.62 2.B - 3.8 
10.0 - 10.79 2.6 - 4.5 


Table 2. Accuracy. 


Resolution 

Ineffectiveness 

At • 0.001 0 
As - 0.1 

0.207 

At «* 0.001 a 
As 10 0.05 

0.202 

At » 6.0005 8 
As • 0.025 

0.201 

At c 0.004 a 
Ax « 0.05 

0.23* 


*Did not reach a steady state — 
t&asa loss 1? per cycle. •- 


In table 3 we ccspare the results for esse I obtained from the two node Is. The full equation 
solution 8ho&?3 a email oscillation in the velocity field which Is not present in the solution using 
the reduced equations. The CPU tiae is on a CSC Cyber 750. The results for the two models sees to 
be in good agreement. 

5.4 Achieving a quasi-steady state 

The results described above were obtained by running the model for 66 cycles. At this point 
there is very little change of average or cajclssura values froa one cycle to the next. For example, 
in case I the EaKimua gas temperature at the cold end io changing at a rate of I0“ 4 K per cycle. 

This is a very email change; however* this rate appears to be virtually constant over the last 10 or 
20 cycles. Therefore, we don’t know how close we are to a steady oscillation. At this rate 1000 
cycles would be required to effect a 10 percent change in this maximum temperature. We need to find 
& oethod to accelerate convergence. 


Table 3. Comparison of Seduced and Full Equations Models* 


Model 

Ineffectiveness 

Temperature 

Swing 

K 

Pressure 

Swing 

MPa 

CPC 

Tice/Cycle 

B 

Reduced 

0.202 

10.0 - 11.2 

2.7 - 4.3 

6.8 

Full 

— ... - 

0.136 

10.0 - XI. 1 

2.8 - 4.3 

23 
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Notation 



SI Units 

P 

- 

pressure 

KPa 

0 

- 

gas density 

Kg/© 3 

Pm 

- 

rntrix density 

kg/n 3 

c 

- 

acoustic velocity 

m/s 

u 

- 

velocity 

ta/s 

G r 

- 

Crueneieen parameter 


h 

- 

heat transfer coefficient 

W/o^K 

*>h 

- 

hydraulic diameter 

Q 

T 

- 

gas temperature 

k 


- 

matrix teaperature 

k 

F 

- 

friction factor 


Cv 

- 

gas heat capacity 

J/kg*l< 

°ta 

- 

matrix heat capacity 

J/kg-X 

<> 

- 

porosity 



- 

gas constant 

J/kg-K 

Y 

- 

ratio of specific heats 


L 

- 

length of regenerator 

Q 

a(t) 

- 

left endpoint (piston position) 

B 

b(t) 

- 

right endpoint 

B 
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